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Inner  shell  atomic  electrons  that  are  tunnel  ionized  in  multi-petawatt  class  laser 
pulses  are  accelerated,  in  vacuum,  to  multi-GeV  energies  in  the  forward  direction. 
In  extreme  fields,  tunnel  ionized  electrons  can  be  brought  to  the  speed  of  light  so 
abruptly,  that  they  stay  in  the  same  phase  of  the  laser  field  throughout  significant 
portions  of  the  confocal  region.  An  analysis  of  the  acceleration  process  is  given,  and 
relativistically  covariant  four-dimensional  numerical  calculations  especially  suited 
for  extreme  fields  are  carried  out.  Radiation  reaction  is  included,  and  the  latest 
relativistic  tunneling  ionization  theories  are  used  to  spawn  the  simulated  electrons. 
An  experimental  configuration  is  suggested,  utilizing  the  10  petawatt  ELI-NP  laser, 
and  plasma  lens  assisted  focusing. 
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I.  INTRODUCTION 


This  report  concerns  a  mechanism  of  electron  acceleration  by  free  space  electromagnetic 
fields,  wherein  the  electrons  are  supplied  by  the  inner  shell  orbitals  of  moderately  heavy 
atoms.  The  process  involves  both  relativistic  tunneling  ionization,  and  classical  charged 
particle  dynamics  in  ultrarelativistic  fields.  With  the  imminent  commissioning  of  the  10 
petawatt  (PW)  beamlines  at  the  ELI-NP  facility  [1],  the  proposed  scheme  promises  to 
advance  the  state  of  the  art  of  free  space  acceleration  of  electrons  by  about  3  orders  of 
magnitude. 

Acceleration  of  electrons  by  electromagnetic  fields  in  free  space  has  been  investigated  at 
length  over  the  years  [2-8].  Free  space  acceleration  is  usefully  framed  in  terms  of  the  Lawson- 
Woodward  (LW)  theorem  [9],  which,  roughly  speaking,  states  that  upon  linearizing  the 
forces,  the  net  energy  gain  is  zero.  Conventional  accelerators  overcome  LW  by  introducing  a 
metallic  structure  comparable  in  dimension  to  the  wavelength  of  the  electromagnetic  field. 
Schemes  such  as  the  inverse  Cherenkov  accelerator  [2]  overcome  LW  by  performing  the 
acceleration  in  a  gaseous  medium  (i.e.,  not  in  free  space).  Several  schemes  rely  on  the 
ponderomotive  force  to  supply  a  nonlinearity  which  overcomes  LW  [3,  5,  7].  In  the  case  of 
Laser  Ionization  and  Ponderomotive  Acceleration  (LIPA)  [5],  an  additional  consideration 
is  that  tunnel-ionized  electrons  are  introduced  into  the  high-field  region  abruptly,  further 
stressing  the  assumptions  of  LW.  In  the  ultrarelativistic  limit,  the  nonlinear  forces  cannot  be 
described  ponderomotively,  because  of  the  possibility  of  a  slow  rate  of  phase  slippage.  This 
is  the  basis  of  the  “Capture  and  Accelerate  Scenario”  (CAS)  [6],  for  which  there  is  recent 
experimental  evidence  [8]. 

This  report  analyzes  and  simulates  LIPA  in  extreme  fields  (xLIPA).  In  the  original  LIPA 
scheme,  electrons  originate  from  a  tenuous  gas  of  moderately  heavy  elements.  LIpon  focusing 
a  high-intensity  laser  pulse  into  the  gas,  electrons  are  tunnel  ionized  and  accelerated.  It  is 
desirable  in  the  LIPA  scheme  for  the  ionization  potential  to  be  matched  to  the  laser  power 
and  focusing  such  that  the  electrons  are  ionized  only  when  they  are  exposed  to  the  peak  laser 
intensity.  Such  electrons  experience  the  highest  ponderomotive  potential,  and  therefore  gain 
the  greatest  energy  (electrons  ionized  earlier  in  the  pulse  may  be  pushed  out  of  the  focus 
before  the  peak  of  the  pulse  arrives).  In  the  case  of  xLIPA,  the  ionization  potential  has 
to  be  matched  in  a  similar  way,  but  the  acceleration  is  no  longer  ponderomotive.  Instead, 
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electrons  stripped  from  atoms  with  certain  favorable  initial  positions  are  brought  abruptly 
to  the  speed  of  light  and  gain  further  energy  as  they  stay  in  nearly  the  same  phase  of  the 
laser  field  for  extended  periods  of  time.  This  is  essentially  the  acceleration  mechanism  of 
CAS,  except  that  the  electrons  are  bound  to  a  parent  ion,  and  hence  nearly  immobile,  until 
the  most  intense  portion  of  the  laser  pulses  passes  by.  In  practical  terms,  xLIPA  accesses 
the  CAS  regime  without  the  need  for  a  source  of  externally  injected  electrons,  which  is  a 
significant  experimental  advantage. 

A  schematic  of  the  proposed  experimental  configuration  is  shown  in  Fig.  1.  In  considering 
a  10  PW  laser  system,  one  confronts  the  practical  issue  that  the  final  focusing  optic  is  likely 
to  be  a  large,  one-of-a-kind  parabolic  mirror,  fixing  the  f-number  for  all  experiments.  To 
increase  the  flexibility  of  the  apparatus,  a  plasma  lens  [10-13]  can  be  introduced.  The  plasma 
lens  can  be  positioned  in  regions  where  the  laser  intensity  is  high,  and  so  does  not  have  to 
be  particularly  large.  In  fact,  in  the  case  of  xLIPA,  the  plasma  lens  itself  can  provide  the 
ions  from  which  electrons  are  stripped  and  accelerated.  In  particular,  the  parameters  can  be 
arranged  so  that  the  beam  waist  occurs  near  the  back  surface  of  the  lens,  where  a  tenuous 
plasma  of  partially  ionized  gas  (say,  argon)  awaits  further  ionization  by  the  focused  laser 
pulse.  For  sufficient  focusing  strength,  the  K-shcll  electrons  are  tunnel  ionized  and  brought 
to  GeV  energies.  Apart  from  flexibility,  the  plasma  lens  offers  three  additional  benefits. 
First,  it  is  likely  that  a  plasma  lens  target  counteracts  the  tendency  of  pre-pulses  to  increase 
the  waist  size  of  the  main  pulse.  Second,  the  plasma  lens  acts  as  an  axially  extended  target 


FIG.  1:  Schematic  xLIPA  Configuration.  The  plasma  lens  is  both  the  final  focusing  element  and 
the  target.  It  serves  to  (i)  provide  flexibility  in  focusing  conditions,  (ii)  decrease  sensitivity  to 
pre-pulses  and  pointing  stability,  and  (iii)  provide  an  optimal  configuration  of  ions  from  which 
xLIPA  electrons  are  extracted. 


2 


with  focusing  properties  throughout,  which  reduces  sensitivity  to  errors  in  the  position  of  the 
vacuum  focus.  Finally,  as  will  be  seen  below,  the  highest  energy  xLIPA  electrons  originate 
off-axis,  and  therefore  a  target  with  an  on-axis  density  minimum  leads  to  a  more  favorable 
final  energy  distribution. 


II.  ANALYSIS  OF  THE  ULTRA-RELATIVISTIC  CASE 


There  is  an  exact  solution  for  the  motion  of  a  charged  particle  in  any  radiation  held  of  the 
form  F/1!y(h/tx/J),  where  k/t  is  the  four  dimensional  wavevector  of  the  radiation,  and  x M  is  the 
spacetime  coordinate.  Without  loss  of  generality,  let  k\  —  k-z  —  0.  Then  uq  —  u^  is  invariant, 
where  u M  is  the  four  dimensional  velocity  of  the  particle1.  Combining  this  with  the  identity 
=  1,  one  can  show  that  in  the  high  energy  limit  u3  3>  {ui,^}.  In  other  words, 
the  particle  moves  predominatly  in  the  “forward”  direction,  i.e.,  parallel  to  the  radiation 
wavevector. 

The  following  analysis  is  based  on  the  expectation  that  in  extreme  laser  holds,  an  inner 
shell  electron  that  is  tunnel- ionized,  is  accelerated  abruptly  to  the  speed  of  light.  Then,  ac¬ 
cording  to  the  foregoing  discussion,  the  direction  of  motion  is  nearly  parallel  to  the  wavevec¬ 
tor  of  the  radiation,  and  the  phase  of  the  particle  in  the  radiation  held  can  be  regarded 
as  constant.  The  primary  constraint  is  that  the  interaction  is  limited  to  regions  where  the 
irradiance  is  high  and  the  phase  velocity  is  close  to  c.  This  corresponds  to  the  two  regions 
just  outside  the  confocal  region.  That  is,  far  from  the  confocal  region  the  irradiance  is  too 
low,  but  inside  the  confocal  region  the  phase  velocity  is  too  high. 

The  exact  equations  of  motion  for  a  particle  in  a  plane  wave,  written  as  matrix  equations, 


are 


dx 

—  =  cu 
ds 

ds 


(1) 

(2) 


Here,  x(s)  is  the  world  line  of  the  particle,  u(s)  is  the  four- velocity,  and  U  =  a(s)uF.  The 
parameter  s  is  the  proper  time,  a(s)  =  qE(s)/mcou,  E(s )  is  the  electric  held,  q  is  the  charge 
of  the  particle,  m  is  the  mass,  and  ui  is  the  frequency  of  the  radiation.  Using  the  coordinate 


1  In  particular,  u  =  (7,  7/3i,  7/?25  7/?3)i  where  c/3?:  are  the  components  of  the  three  dimensional  velocity,  and 
7=  (1  -P2)~1/2 
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system  described  above,  the  matrix  F  is 

/  0  1  0  0 

10  0-1 

F  = 

0  0  0  0 

\  0  1  0  0 

As  long  as  the  particle  stays  in  phase,  hi  can  be  regarded  as  constant,  and  the  solution  of 
the  velocity  equation  is 

u(s)  =  ensw(0)  (4) 

The  matrix  exponential  ens  is  easily  calculated  due  to  the  nilpotency  of  F .  i.e.,  F3  =  0.  The 
result  is 

1  +  a2  j2  a  0  — a2/ 2 

cr  10  — a 

0  0  1  0 

<72/2  a  0  1  —  a2 /2 

where  a(s)  =  aus.  ft  is  easily  verihed  that  A  is  a  Lorentz  transformation,  i.e.,  ATgA  =  g, 
where  T  indicates  the  transpose,  and  g  =  diag(l,  —1,  —1,  —1).  Of  particular  interest  is  the 
initial  condition  w(0)  =  (1,0,0,0)T,  which  according  to  most  theories  holds  for  an  electron 
at  the  moment  of  ionization,  at  least  when  the  atomic  number  satisfies  Z  <C  137.  In  this 
case, 

/  1  +  a2 12  \ 

«(s)  =  (6) 

0 

V  cr2 /2 

Note  that  the  invariant  Uq  —  is  maintained.  Moreover,  when  a  =  aus  3>  1  the  momentum 
is  predominantly  in  the  forward  direction,  i.e.,  M3  U\.  Assuming  a>  1,  this  requires  that 
us  be  at  least  of  order  unity,  i.e.,  the  time  elapsed  according  to  a  clock  moving  with  the 
particle  should  read  at  least  one  laser  period,  as  measured  by  a  lab  frame  clock.  This  does 
not  necessarily  violate  the  assumption  that  the  particle  should  stay  in  phase,  since  the  two 
clocks  may  keep  very  different  time. 

Since  the  Rayleigh  length  characterizes  the  interaction  length,  which  is  expressed  in  terms 
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of  x3,  it  is  useful  to  put  s  in  terms  of  x3.  These  are  related  by 


x3  =  c  /  u3(s)ds  =  -a2u)2s3 

J  6 

Defining  k  =  co/c  gives 

(6a;3)1/3 
S  c(ak)2/3 

Then  the  energy  of  the  particle  is 

u0(x3)  =  1  +  i  (6akx3)2/3 


(7) 

(8) 

(9) 


Here,  it  is  understood  that  the  spacetime  origin  is  chosen  so  that  x3  =  0  at  the  moment 
of  ionization.  Finally,  the  highest  possible  energy  is  estimated  by  substituting  the  Rayleigh 
length  for  x3,  i.e.,  x3  — >  vr^/A,  where  r0  is  the  radius  of  the  beam  waist,  and  A  =  2n/k. 
This  gives 

''-Umax  =  1  +  ^  (x)4/3  (12fl)2/3  (10) 

or 

«0,max~12(y)4/3a2/3  (11) 

As  an  example,  a  10  PW  laser  pulse,  with  A  =  0.8  /ini,  focused  to  ro  =  5  /mi,  gives  a  ~  100, 
and  Mo,max'nrc2  ~  1.5  GeV.  The  question  of  how  to  match  the  ionization  potential  of  the  ion 
to  the  laser  parameters,  so  that  the  electron  becomes  free  at  the  optimal  time,  is  addressed 
below. 


III.  NUMERICAL  SIMULATIONS 

The  primary  observable  in  an  xLIPA  configuration  is  the  electron  momentum  distribution. 
Under  the  assumption  that  the  electrons  are  drawn  from  a  gas  whose  density  is  low  enough 
so  that  plasma  effects  are  negligible,  the  momentum  distribution  can  be  obtained  by  single 
particle  tracking  methods.  In  the  case  of  extreme  fields,  obtaining  an  accurate  numerical 
solution  to  the  equation  of  motion  is  non-trivial.  The  methods  used  in  this  work  are  detailed 
in  section  IV.  In  addition  to  integrating  the  equation  of  motion,  the  calculation  has  to 
account  for  tunneling  ionization.  For  this  purpose,  an  analytical  rate  law  is  needed.  In  this 
work,  the  Coulomb-corrected,  dressed  strong  field  approximation  (SFA)  of  Klaiber  et  al. 
[14]  is  used.  The  particle  tracking  code  used  here  also  allows  the  exact  Landau  and  Lifshitz 
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radiation  reaction  formula  to  be  incorporated  into  the  integration.  It  is  found  that  at  the  10 
PW  level  this  effect  is  very  minor.  An  example  based  on  more  speculative  laser  parameters 
is  given  where  the  effect  is  noticeable. 

A.  Choice  of  Ionic  Species  and  Ionization  Thresholds 

As  is  well  known,  the  rate  of  tunneling  ionization  as  a  function  of  the  electric  held  is 
exponential,  so  that  the  concept  of  a  threshold  held  is  applicable.  The  threshold  for  a 
given  ionization  potential  determines  where  in  the  laser  focus  a  free  electron  hrst  appears. 
The  ionization  potential  that  leads  to  the  highest  energy  electrons  has  to  be  calculated  by 
performing  a  series  of  numerical  experiments  for  a  given  laser  power  and  focusing  configu¬ 
ration.  One  expects  that  there  is  a  finite  optimum,  for  a  threshold  that  is  too  high  restricts 
the  starting  coordinate  to  a  small  region  near  the  focus,  while  one  that  is  too  low  restricts 
the  starting  coordinate  to  a  region  far  from  the  focus,  where  the  fields  are  too  small  to 
accelerate  electrons  in  the  forward  direction  before  ejecting  them  radially. 

One  useful  estimate  for  the  threshold  held  is  the  one  due  to  Augst  [15],  which  is  based 
on  a  simple  barrier  suppression  picture,  and  is  in  reasonable  agreement  with  experimental 
data.  An  alternative  threshold  based  on  the  relativistic,  Coulomb-corrected  SFA  (or  any 
rate  law)  can  be  derived  as  follows.  Let  the  rate  be  given  as  some  function  w(Ui,  E0),  where 
Ui  is  the  ionization  potential,  and  E0  is  the  peak  electric  held  of  the  applied  laser  pulse. 
The  threshold  held  is  that  value  of  E0  which  satishes  w(Ui,E0 )  =  u/27t,  where  u  is  the 
laser  frequency.  Note  that  since  the  dependence  of  the  left  hand  side  on  Eq  is  generally 
exponential,  the  result  only  depends  logarithmically  on  the  choice  of  time  scale  appearing 
on  the  right  hand  side. 

Ionization  thresholds  for  various  ions  are  displayed  in  Table  I  based  on  the  dressed 
Coulomb-corrected  SFA  [14] .  The  charge  states  displayed  are  selected  based  on  their  status 
as  the  highest  charge  states  for  a  given  set  of  non-empty  shells,  e.g.,  Ar'+  is  the  highest  ar¬ 
gon  charge  state  that  still  has  an  electron  in  the  M  shell.  When  the  atom  is  first  exposed  to 
the  laser  radiation,  electrons  are  expected  to  be  stripped  sequentially,  in  order  of  increasing 
ionization  potential.  Therefore  the  potentials  used  at  each  stage  are  assumed  to  be  those 
appropriate  for  an  ion  in  the  lowest  energy  state.  It  should  be  noted  that  the  SFA  rate 
law  is  strictly  valid  only  for  hydrogen-like  ions.  In  addition  to  the  intensity  threshold,  the 
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TABLE  I:  Tunneling  Ionization  Threshold  for  Ions  of  Interest 


Ion 

Potential  (eV) 

Occ.  Shells 

SFA  Threshold0  (W/crn2) 

4hresh  (A  =  800  nm) 

Ar7+ 

143 

KLM 

1.4  x  1017 

0.26 

Ar15+ 

918 

KL 

2.7  x  1019 

3.57 

Ar17+ 

4426 

K 

2.4  x  1021 

33.6 

Ti21+ 

6628 

K 

7.6  x  1021 

59.8 

Fe25+ 

9278 

K 

2.0  x  1022 

96.9 

J£r35+ 

17948 

K 

1.3  x  1023 

250 

Xe53+ 

41347 

K 

1.5  x  1024 

834 

Au78+ 

93459 

K 

1.6  x  1025 

2740 

U91+ 

132280 

K 

4.5  x  1025 

4570 

“The  static  field  rate  is  used  to  estimate  the  threshold.  Cycle  averaging  increases  the  threshold  by  «  25%. 

corresponding  peak  normalized  vector  potential,  ao,  is  displayed.  The  parameter  a0  is  of 
fundamental  importance  because  the  interaction  becomes  ultra-relativistic  when  ao  3>  1. 

The  laser  pulse  format  assumed  throughout  section  III  is  a  waist  radius  (1/e  of  the  field) 
of  5  jum,  a  pulse  length  (1/e  of  the  field)  of  20  fs,  a  wavelength  of  0.8  jum,  and  a  field 
configuration  consistent  with  the  lowest  order  Hermite-Gaussian  mode.  Longitudinal  field 
components  are  chosen  to  be  consistent  with  the  Coulomb  gauge  (V  ■  A  =  0),  including 
a  correction  to  account  for  finite  pulse  length.  Two  laser  powers  of  primary  interest  are 
considered:  10  PW  giving  a0  =  100,  and  25  exawatts  (EW)  giving  o0  =  5000.  The  lat¬ 
ter  stupendous  laser  power  is  chosen  because  it  gives  a  remarkable  signature  of  radiation 
reaction.  The  former  is  chosen  because  it  corresponds  to  the  near-term  ELI-NP  parameters. 

The  effect  of  the  ionization  potential  on  the  xLIPA  distribution  for  the  10  PW  case  is 
shown  in  Fig.  2.  The  ionization  potentials  used  correspond  to  5  different  ionic  species,  Ar'  +  , 
Ar15+,  Ar17+,  Ti21+,  and  Fe25+.  The  latter  three  ions  are  all  hydrogen-like.  Representative 
particle  orbits  from  the  5  cases  are  shown  in  Fig.  2(a),  in  the  plane  of  energy  and  the 
longitudinal  coordinate,  X3.  The  coordinate  is  raised  to  the  1/3  power  in  order  to  achieve  the 
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(a) 


(b) 


Ionization  Potential  (eV) 


FIG.  2:  Effect  of  ionization  potential  on  electron  acceleration  for  ao  =  100.  Representative  orbits 
for  ionization  of  5  different  ionic  species  are  shown  in  (a),  while  the  highest  energy  selected  from 
an  initial  distribution  of  106  particles  is  shown  in  (b)  as  a  function  of  ionization  potential,  for  the 
same  5  species. 

effect  of  a  logarithmic  plot,  while  still  allowing  for  a  sign  and  a  zero-crossing2.  The  lowest  two 
charge  states,  with  the  smallest  ionization  potentials,  are  seen  to  experience  ponderomotive 
acceleration,  evidenced  by  the  fast  oscillations  in  energy  superimposed  on  the  overall  energy 
gain.  The  net  energy  gain  is  also  much  smaller  than  that  of  the  higher  charge  states,  and 
the  acceleration  is  terminated  due  to  radial  expulsion  from  the  confocal  region  (not  shown). 
In  contrast  the  three  higher  charge  state  species  exhibit  a  large  energy  gain  within  a  single 
optical  cycle,  as  evidenced  by  a  steep  initial  slope  devoid  of  any  oscillatory  features. 

The  anticipated  peak  in  the  maximum  accelerated  energy  vs.  ionization  potential  is 
illustrated  in  Fig.  2(b).  In  order  to  produce  the  data,  106  particles  are  loaded  into  a  uniform 
distribution  in  a  region  sufficiently  large  to  encompass  the  ionized  volume.  The  particle 
with  the  highest  energy,  evaluated  far  from  the  laser  focus,  is  used  to  define  the  maximum 
accelerated  energy.  In  the  case  a0  =  100,  with  the  assumed  focusing  conditions,  hydrogen¬ 
like  argon  turns  out  to  be  an  optimal  species. 

2  Strictly  this  requires  the  ad-hoc  relationship  (—  |cc|)1//3  =  — |a:|1//3 
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FIG.  3:  Final  distribution  of  electrons  drawn  from  the  K-shell  of  argon  with  10  PW  of  laser  power. 
The  momentum  space  distribution  is  shown  in  (a,b,c),  and  the  angular-spectral  distribution  is 
shown  in  (d).  7/3*  is  momentum  normalized  to  me,  with  laser  polarization  in  the  1-direction  and 
central  wavevector  in  the  3-direction.  The  color  scale  is  logarithmic.  The  dashed  line  in  (d)  is  the 
curve  on  which  all  particles  would  lie  in  the  plane- wave  limit. 

B.  Electron  Distributions  at  muliti-Petawatt  Scale 

Based  on  the  forgoing  discussion,  the  ions  Ar16+  and  Ar1,+  (both  with  similar  ionization 
potentials)  are  of  interest  in  terms  of  possible  near-term  ELI-NP  experiments  utilizing  the 
10  PW  beamlines.  In  order  to  broadly  characterize  the  outgoing  electron  distribution, 
it  is  sufficient  to  uniformly  load  the  ions  of  interest  into  a  box  encompassing  the  focal 
volume.  Then,  by  examining  the  correlation  between  final  energy  and  initial  position  (or 
other  classical  5-matrix  projections),  more  favorable  distributions  can  be  identified.  In  an 
experiment,  there  is  no  need  to  pre-form  Ar16+  or  Ar17+  in  the  focal  volume.  So  long  as 
any  lower  charge  state  of  argon  is  present,  and  provided  the  main  10  PW  pulse  comes  to 
a  suitable  focus,  all  higher  charge  states  will  be  sequentially  produced  as  the  applied  held 
increases.  The  use  of  a  plasma  lens  to  assist  the  focusing  is  addressed  in  section  V. 
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FIG.  4:  Classical  S-matrix  projections  for  xLIPA  electrons  in  the  10  PW  case.  The  vertical  axes  are 
the  final  particle  energies  and  the  horizontal  axes  are  the  initial  positions  x*,  with  laser  polarization 
in  the  1-direction  and  central  wavevector  in  the  3-direction.  High  energies  are  encouraged  by 
increasing  the  ion  density  near  the  two  positions  (xi,  X2,  x3)  ~  (±7,  0,  —200)  /irn. 

The  electron  distribution  from  a  uniform  background  of  Ar17+  is  shown  in  Fig.  3.  The 
momentum  distribution  is  given  in  the  three  possible  two-dimensional  momentum  planes. 
A  strong  preference  for  the  polarization  plane  is  indicated  in  Fig.  3(c),  given  that  the  color 
scale  is  logarithmic.  The  angle-energy  plane  is  a  promising  experimental  observable.  A 
distinct  population  of  high  energy  electrons  can  be  seen  around  1.5  GeV,  at  an  angle  of  a 
few  milliradians  with  respect  to  the  laser  propagation  axis.  In  a  separate  simulation  with 
the  longitudinal  field  components  artificially  suppressed,  the  energy  gain  was  about  half  that 
obtained  with  the  longitudinal  fields  included. 

It  is  possible  to  optimize  the  electron  energy  distribution  by  controlling  the  initial  spatial 
distribution  of  the  ions.  In  order  to  determine  how  to  do  this,  a  classical  S’-matrix  is 
computed  by  correlating  the  initial  and  final  states  of  a  large  number  of  particles.  Projections 
of  this  S-matrix  are  shown  in  Fig.  4.  In  interpreting  the  figure,  the  limitations  of  two- 
dimensional  plotting  have  to  be  taken  into  account.  The  mapping  between  the  initial  spatial 
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point  (xi,X2,xs)  and  the  final  energy  is  much  closer  to  one-one  than  can  be  gleaned  by 
inspection  of  any  of  the  three  projections.  The  only  source  of  dispersion  in  the  spectrum 
of  particles  originating  from  a  single  spatial  point  is  the  variability  in  the  phase  at  the 
moment  of  ionization,  which  in  turn  is  clue  to  the  fundamental  statistical  nature  of  quantum 
mechanics. 

The  optimal  distribution  of  ions,  based  on  Fig.  4,  is  the  pair  of  points  (aq,  X2,  X3)  ~ 
(±7, 0,  —200)  /am.  An  approximation  of  this  would  be  a  ring  with  diameter  14  /am,  centered 
on  the  £3  axis,  and  upstream  of  the  laser  focus  by  200  /am.  A  further  approximation  of  this 
would  be  a  plasma  lens,  which  has  a  natural  density  minimum  on-axis,  situated  with  the 
laser  focus  slightly  beyond  the  lens  exit  plane.  Thus,  as  mentioned  above,  a  plasma  lens  has 
the  fortuitous  property  of  weighting  the  electron  spectrum  toward  higher  energies. 

C.  Electron  Distributions  at  multi-Exawatt  Scale 

In  this  section,  the  focusing  conditions  are  the  same  as  in  the  10  PW  case  discussed 
above,  but  the  laser  power  is  increased  to  25  EW.  This  gives  ao  =  5000,  which  suggests, 
based  on  Table  I  and  Fig.  2,  that  the  K-shcll  of  gold  is  a  suitable  source  of  electrons  in 
this  case.  Although  present  day  laser  technology  is  far  from  producing  25  EW  of  power, 
it  is  interesting  to  investigate  how  the  xLIPA  mechanism  scales,  and  whether  unambiguous 
signatures  of  radiation  reaction  are  obtained. 

Consider  first  the  case  with  radiation  reaction  neglected.  The  final  electron  distributions 
at  25  EW,  analagous  to  the  10  PW  distributions  of  Fig.  3,  are  shown  in  Fig.  5.  The  most 
obvious  difference  between  the  two  cases  is  that  the  highest  energy  particles  are  on-axis,  and 
the  highest  energy  is  about  70  GeV  rather  than  2  GeV.  Linear  scaling  with  ao  would  give 
100  GeV,  while  the  scaling  of  Eq.  (11)  would  give  30  GeV. 

The  S-matrix  projections  at  25  EW,  analagous  to  the  10  PW  projections  of  Fig.  4,  are 
shown  in  Fig.  6.  While  there  are  some  similarities,  a  notable  difference  is  that  some  high 
energy  particles  originate  on-axis.  This  may  be  an  indication  that  the  simple  plane  wave 
analysis  of  section  II  holds  more  closely  in  the  25  EW  case. 
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FIG.  5:  Final  distribution  of  electrons  drawn  from  the  K-shell  of  gold  with  25  EW  of  laser  power, 
neglecting  radiation  reaction.  The  momentum  space  distribution  is  shown  in  (a,b,c),  and  the 
angular-spectral  distribution  is  shown  in  (d).  7 0,  is  momentum  normalized  to  me,  with  laser 
polarization  in  the  1-direction  and  central  wavevector  in  the  3-direction.  The  color  scale  is  log¬ 
arithmic.  The  dashed  line  in  (d)  is  the  curve  on  which  all  particles  would  lie  in  the  plane-wave 
limit. 


D.  Radiation  Reaction 

Radiation  reaction  (RR)  is  one  of  the  few  areas  of  classical  electrodynamics  that  is  still 
lacking  in  experimental  confirmation,  mainly  due  to  the  lack  of  experimental  facilities  ca¬ 
pable  of  accessing  the  radiation  dominated  regime.  This  issue  has  attracted  significant 
attention  in  recent  years  [16-21].  The  xLIPA  process  is  a  candidate  for  experimental  obser¬ 
vation  of  RR.  In  particular,  if  one  can  show  that  the  angle-energy  distribution  is  affected  in 
a  clear  manner  by  RR,  this  might  lead  to  a  powerful  confirmation,  or  denial,  of  the  valid¬ 
ity  of  the  existing  theories.  The  particle-tracking  code  used  in  this  work  has  the  option  of 
including  the  Landau  and  Lifshitz  RR  formula  [22],  It  reproduces  very  closely  the  results 
reported  in  Ref.  [21]. 
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FIG.  6:  Classical  S'-matrix  projections  for  xLIPA  electrons  in  the  25  EW  case,  neglecting  radiation 
reaction.  The  vertical  axes  are  the  final  particle  energies  and  the  horizontal  axes  are  the  initial 
positions  Xi,  with  laser  polarization  in  the  1-direction  and  central  wavevector  in  the  3-direction. 
The  vertical  feature  at  x±  =  0  in  (a)  is  an  artifact  of  the  initial  particle  loading. 

Before  addressing  the  effect  of  RR  on  xLIPA,  it  is  useful  to  consider  the  effect  on  the 
well-known  solution  for  the  motion  of  an  electron  in  a  plane-wave.  This  is  illustrated  in 
Fig.  7,  for  the  case  a0  =  1000.  Figs.  7(a,c,e)  show  the  orbits  neglecting  RR,  as  computed 
by  the  particle-tracking  code.  The  numerical  error  can  be  gauged  by  tracking  the  invariant, 
T  =  Mo  —  U3,  which  for  the  chosen  initial  conditions  should  be  unity  at  all  times.  The 
observed  error  is  less  than  1%.  Figs.  7(b,d,f)  show  the  same  orbits,  except  that  RR  is 
accounted  for  using  the  exact  Landau  and  Lifshitz  theory.  The  primary  difference  is  that 
the  formerly  invariant  T  decreases  at  each  turning  point  in  the  electron  motion.  This  is 
consistent  with  expectations,  since  the  turning  points  are  where  the  acceleration,  and  hence 
radiated  power,  are  greatest.  Interestingly,  a  close  look  at  Fig.  7(d)  shows  that  M3  reaches  a 
higher  value  during  the  second  half-cycle  of  the  motion  compared  with  the  first.  This  implies 
that  the  energy  also  increases,  since  Wo  =  M3  +  T,  and  T  changes  only  slightly.  Naturally, 
energy  need  not  be  conserved  in  the  presence  of  an  external  held. 
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FIG.  7:  Electron  orbits  in  a  plane  wave  without  the  reaction  force  (a,c,e)  and  with  the  exact 
Landau  and  Lifshitz  reaction  force  (b,d,f).  Here  T  =  uq  —  u%,  which  is  expected  to  be  invariant  in 
the  absence  of  radiation  reaction. 

Fig.  8  shows  a  set  of  multi- EW  xLIPA  distributions  with  RR  included.  These  should  be 
compared  with  Fig.  5,  which  is  an  identical  case  with  RR  neglected.  The  overall  effect  of 
RR  is  to  narrow  the  distributions,  and  to  make  the  cutoffs  more  abrupt.  Considering  the  log 
scale,  Fig.  8(a)  and  (c)  indicate  that  in  the  presence  of  RR,  there  are  two  fairly  well  defined 
beamlets  with  m2  ~  ±1000  and  u 3  ~  20000.  In  the  absence  of  RR  these  features  become 
much  more  spread  out  in  momentum  space. 


IV.  NUMERICAL  METHODS 

The  numerical  model  employed  in  this  work  incorporates  a  number  of  advances,  includ¬ 
ing  a  covariant  particle  pusher  that  respects  Lorentz  invariance  to  machine  precision,  and 
elegantly  incorporates  radiation  reaction  and  automatic  time  step  adjustment.  The  ioniza¬ 
tion  model  uses  the  latest  relativistic  theories,  and  confines  the  use  of  the  random  number 
generator  to  the  initial  conditions.  Finally,  the  implementation  takes  advantage  of  hardware 
acceleration  by  means  of  general  purpose  graphical  processing  units. 
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FIG.  8:  Final  distribution  of  electrons  drawn  from  the  K-shell  of  gold  with  25  EW  of  laser  power,  in¬ 
cluding  radiation  reaction.  The  momentum  space  distribution  is  shown  in  (a,b,c),  and  the  angular- 
spectral  distribution  is  shown  in  (d).  7/%  is  momentum  normalized  to  me,  with  laser  polarization 
in  the  1-direction  and  central  wavevector  in  the  3-direction.  The  color  scale  is  logarithmic.  The 
dashed  line  in  (d)  is  the  curve  on  which  all  particles  would  lie  in  the  plane- wave  limit,  in  the 
absence  of  RR. 


A.  Covariant  Particle  Pusher 

The  particle  pusher  calculates  the  world  lines  of  particles  introduced  into  the  external 
field.  Countless  authors  have  implemented  schemes  for  this  purpose.  One  of  the  most 
commonly  used  schemes  is  the  one  introduced  by  Boris  [23].  I11  carrying  out  simulations  of 
laser-particle  interactions,  one  most  often  chooses  a  time  step  that  is  some  fixed  fraction  of 
the  laser  period.  In  the  present  case,  this  is  not  satisfactory  because  in  the  ultrarelativistic 
limit,  there  are  turning  points  where  the  characteristic  time  scale  is  a  very  small  fraction 
of  the  laser  period.  At  the  same  time,  one  must  propagate  each  particle  to  a  point  well 
outside  the  confocal  region.  These  two  requirements  introduce  a  large  scale  separation, 
which  demands  some  form  of  automatic  time  step  adjustment.  It  turns  out  that  this  issue 
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can  be  resolved  elegantly  by  designing  the  pusher  to  operate  in  covariant  fashion,  i.e.,  to 
push  particles  in  proper  time.  In  addition,  a  covariant  particle  pusher  allows  the  exact 
Landau  and  Lifshitz  radiation  reaction  formula  to  be  incorporated  in  a  simple  way. 

As  in  the  Boris  pusher,  we  choose  to  leapfrog  momentum  and  position,  except  that 
these  are  now  four-vectors,  with  the  independent  variable  being  proper  time.  The  position 
equation  can  be  updated  trivially  using 

x(s  +  As)  =  x(s)  +  cu(s  +  As/2)As  (12) 


where  x  and  u  are  four-element  column  vectors.  The  matrix  equation  for  the  momentum  is 
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Here,  the  fields  E  and  B  are  considered  functions  of  s.  If  they  are  prescribed  as  functions 
of  the  spacetime  coordinates,  as  is  typical,  then  one  uses,  e.g.,  the  functional  composition 
E(s)  =  E(  x)  o  x(s),  where  x(s)  is  the  solution  of  (12).  Clearly,  any  valid  must  be  the 
generator  of  a  Lorentz  transformation,  since  uTu  =  1  is  an  identity. 

Assume  that  the  time  step  As  is  chosen  to  be  small  enough  so  that  is  nearly  constant. 
Then 

u(s  +  As)  =  A(s  +  As/2,  As)u(s)  (15) 


where  A(s,  As)  =  en^As.  The  matrix  exponential  of  H  is  unwieldy,  but  can  be  greatly 
simplified  by  performing  a  decomposition  between  electric  and  magnetic  fields.  In  particular, 
using  the  Campbell-Baker-Hausdorf  expansion  for  the  exponential  of  a  sum,  one  obtains  the 
preliminary  result 

A  =  AEABAExB  +  0(As3)  (16) 


where  AE  is  a  boost,  AB  is  a  rotation,  and  AExB  is  a  boost  that  corrects  for  the  fact  that 
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boosts  and  rotations  do  not  commute.  More  specifically, 
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where  y  =  1  —  cos 9,  ij)  =  sin0,  6  =  q\T5\As/mc,  and  b  =  B/|B|.  The  AExB  boost  has 
the  same  form  as  AE,  with  (3  =  q2 |E  x  B|As2/m2c2  and  e  =  E  x  B/|E  x  B|.  However,  it 
turns  out  that  the  explicit  application  of  the  correction  AExB  is  not  the  most  efficient  way 
to  achieve  convergence.  Instead,  one  can  split  the  step  in  a  manner  similar  to  the  Boris 
scheme  to  obtain 


A(s,  As)  ~  Ae(s,  As/2)Ae(s,  As)Ae(s,  As/2)  (19) 


Based  on  a  range  of  numerical  experiments,  this  is  significantly  more  accurate  than  (16). 

In  summary,  updating  u  is  a  matter  of  carrying  out  three  explicitly  given  linear  trans¬ 
formations.  The  update  is  accurate  to  order  As2,  and  respects  Lorentz  invariance,  i.e., 
ATgA  =  g,  to  machine  precision.  One  caveat  is  that  the  unit  vectors  are  ill-defined  in  a 
field- free  region.  In  practice,  this  is  easily  remedied,  e.g.,  by  superposition  of  a  minuscule 
uniform  field  with  the  field  of  interest. 


B.  Timestep  Adjustment,  Accuracy  and  Performance 


A  useful  feature  of  a  covariant  particle  pusher  is  that  a  constant  time  step  (in  proper 
time)  leads  to  a  constant  phase  step  for  any  particle  in  a  plane  wave.  To  see  this,  consider 
the  total  derivative  of  the  phase 


dip  d{k^xh  dx M 


(20) 
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According  to  the  exact  solution  of  the  equations  of  motion,  kn'u/1  is  invariant.  Defining 
T  =  k^u^/cu,  and  demanding  the  that  phase  step  A <p  <C  2tx,  one  obtains  the  necessary 
accuracy  condition 

As  «  %  (21) 

U)  1 

which  involves  only  invariant  quantities.  For  a  particle  initially  at  rest,  such  as  a  tunnel 
ionized  electron,  one  has  T  =  1,  so  that  the  proper  time  step  appropriate  for  a  relativistic 
particle  is  the  same  as  the  lab  frame  time  step  appropriate  for  a  non-relativistic  particle. 
The  reason  for  this  is  that  the  lab  frame  time  step  is  longer  than  the  proper  time  step  by  a 
factor  of  7,  which  is  just  the  right  factor  to  keep  the  phase  step  constant.  It  should  be  noted 
that  in  cases  where  a  laser  pulse  collides  with  an  electron  beam  in  a  counter-propagating 
geometry,  T  can  be  large,  requiring  As  to  be  correspondingly  reduced. 

The  condition  (21)  alone  is  not  sufficient  to  guarantee  accuracy,  because  of  the  Campbell- 
Baker-Hausdorf  expansion  of  the  matrix  exponential.  In  order  to  estimate  the  time  step 
needed  to  justify  truncation  of  the  expansion,  one  may  demand  that  A ExB  should  be  nearer 
the  identity  matrix  than  either  A^  or  A#.  Any  of  these  matrices  approach  the  identity 
as  the  angle  (whether  appearing  in  the  argument  of  a  hyperbolic  or  ordinary  trigonometric 
function)  vanishes.  Therefore,  demanding  that  the  angle  appearing  in  A  exb  should  be  much 
less  than  that  appearing  in  A^  or  A#,  provides  the  required  accuracy  condition.  The  overall 
result  is 
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As  in  the  case  of  advancing  the  particle  momentum,  one  has  an  undefined  floating  point 
operation  in  field-free  regions.  The  remedy  is  again  to  add  a  minuscule  uniform  field 

In  practice  it  is  convenient  to  have  an  expression  for  the  timestep,  As,  not  involving  con¬ 
ditionals,  and  depending  on  a  small  number  of  dimensionless  free  parameters  that  quantify 
the  desired  accuracy.  A  suitable  expression  that  is  used  in  this  work  is 

-1/2 


As  = 


»■(£) 


(E2  +  B2) 


(23) 


The  dimensionless  free  parameters  are  D,  and  1Z,  which  can  be  thought  of  as  frequency 
multipliers.  The  multiplier  D  corresponds  to  the  number  of  steps  taken  during  one  period 
of  the  radiation  held.  The  multiplier  1Z  accounts  for  the  requirement  due  to  non-commuting 
boost  and  rotation  operators. 
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The  effect  of  the  parameter  f2  is  well  known  in  connection  with  linear  interactions,  and 
does  not  require  further  exposition.  The  conservative  value  17  =  267  is  used  throughout 
this  work.  Fig.  9  illustrates  the  accuracy  of  the  pusher  in  the  case  of  a  plane  wave  field 
with  normalized  vector  potential  a(t,z )  =  aosina;(z  —  t).  In  the  given  scenario,  T  =  1  is 
invariant,  so  that  a  suitable  error  measure  is  |1  —  T|,  which  should  vanish  at  all  time  levels. 
The  maximum  value  of  the  error  measure  during  one  period  of  the  motion  is  shown  as  a 
function  of  77,  with  a0  =  100,  in  Fig.  9(a).  Two-digit  accuracy  is  only  obtained  for  77  >  100. 
Using  a  constant  77  =  1000  and  varying  a0  gives  Fig.  9(b).  In  order  to  keep  the  accuracy 
fixed,  one  has  to  keep  77/ao  approximately  fixed.  Available  computation  time  sometimes 
dictates  that  higher  accuracy  is  obtained  for  smaller  a o-  For  the  simulations  discussed  in 
section  III,  we  used  77/ao  =  10  with  a0  =  100,  and  77/ao  =  1  with  a0  =  5000. 

The  covariant  pusher  described  above  is  implemented  as  an  OpenCL  kernel.  The  perfor¬ 
mance  under  various  conditions,  running  on  an  AMD  D700  GPGPU,  is  illustrated  in  Fig.  10. 
The  various  points  on  the  plot  vary  the  number  of  particles  involved  in  the  calculation  (hor¬ 
izontal  axis),  the  floating  point  precision,  and  the  number  of  particle  advances  per  OpenCL 
kernel  invocation.  The  latter  parameter  can  be  important  because  of  the  overhead  involved 
in  invoking  the  kernel,  and  also  because  of  the  cost  of  moving  data  in  and  out  of  cache  at 
the  beginning  and  end  of  the  kernel  invocation,  respectively.  In  the  best  case,  about  1.5  bil¬ 
lion  particles  per  second  can  be  advanced,  including  field  evaluations.  The  double  precision 


(a) 


IE-01 

IE-02 

IE-03 

IE-04 

IE-05 

IE-06 

1E+01  1E+02  1E+03  1E+04 

Normalized  Vector  Potential 


FIG.  9:  Accuracy  of  covariant  pusher,  as  measured  by  the  maximum  value  of  |1  —  T|,  during  one 
period  of  the  motion  of  the  electron  in  a  plane  wave.  The  accuracy  is  displayed  as  (a)  a  function 
of  77,  given  ao  =  100,  and  (b)  as  a  function  of  a o,  given  77  =  1000.  In  all  cases  12  =  267. 
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performance  is  about  1/3  of  this.  The  standard  Boris  pusher  gives  comparable  performance. 


C.  Radiation  Reaction 


A  major  advantage  of  expressing  the  particle  pusher  in  covariant  fashion  is  that  the 
Lorentz-Abraham-Dirac  (LAD)  formula  for  the  radiation  reaction  force  takes  the  simple 
form 


R  = 


2 q2  ( d2u  T  d2u 
3mc  \d^~UU  9lfs 2 


(24) 


The  Landau  and  Lifshitz  (LL)  formula  is  derived  by  substituting  for  d2u/ds 2  the  value  ob¬ 
tained  in  the  absence  of  radiation  reaction.  In  three  dimensional  notation,  the  LL  formula 


is  extremely  unwieldy,  and  even  in  four  dimensional  form,  it  appears  to  require  expensive 
evaulations  of  all  spacetime  derivatives  of  the  held  tensor.  When  the  covariant  pusher  de¬ 
scribed  above  is  used,  a  simple  and  elegant  alternative  becomes  readily  available.  Namely,  by 
splitting  each  step  into  two  half-steps,  d2u/ds2  can  be  evaluated  by  direct  finite  differencing. 
That  is,  during  each  step  generate 
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FIG.  10:  Performance  of  covariant  pusher,  in  terms  of  millions  of  particles  advanced  per  second, 
including  field  evaluations.  In  the  legend,  SP  means  single  precision,  DP  means  double  precision, 
and  s  is  the  number  of  cycles  advanced  for  each  OpenCL  kernel  invocation. 
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This  requires  field  evaluations  at  only  two  spacetime  points.  Now,  in  the  absence  of  radiation 
reaction,  u ^  =  u(s),  u ^  =  u(s  +  As/2),  and  =  u(s  +  As).  The  LL  reaction  force  is 
therefore  given  by  the  matrix  equation 


where 


(S2u  —  uuTgS2u) 

3mc 

(27) 

u(2)  _  2w(l)  +  M(0) 

(28) 

As2/4 

is  the  finite  difference  form  of  d2u/ds2.  The  updated  four- velocity,  including  the  LL  reaction 
force,  is 

u(s  +  As)  =  +  RAs  (29) 

Here,  a  simple  Euler  advance  is  justified  by  the  expectation  that  the  reaction  force  is  small. 


D.  Ionization  Algorithm 


In  the  context  of  a  classical  particle  tracking  calculation,  an  ionization  algorithm  amounts 
to  devising  a  rule  for  spawning  a  particle  in  the  midst  of  the  interaction.  If  the  ion  motion 
is  negligible  on  the  time  scale  of  the  laser  pulse,  an  equivalent  view  is  that  the  ionization 
algorithm  provides  a  rule  for  abruptly  changing  the  charge  of  a  particle  in  the  midst  of  the 
interaction.  In  particular,  the  electron  charge  is  changed  from  q  =  0  to  q  =  —  e  during  one 
time  step. 

In  order  to  avoid  calling  a  random  number  generator  while  integrating  the  equations  of 
motion  (as  is  often  done),  we  associate  with  each  particle  a  constant  parameter  H,  and  an 
evolving  parameter  77(f),  which  encode  all  the  statistical  information  about  the  ionization 
process.  The  condition  for  a  particle  to  be  ionized  is  77(f)  >  H .  The  parameter  H  appears 
in  the  initial  distribution  function 


9(x,P,H)  =  f(x,p)e  H 


where  f(x,p)  is  the  usual  phase  space  distribution.  Let 


(30) 


(31) 


where  7 77(f)  is  the  ionization  rate  evaluated  for  a  given  particle,  which  has  to  be  computed 
using  one  of  several  available  tunneling  theories.  These  usually  have  the  form 

777(f)  =  CiE(t)°2  exp[— Cs/E(t)\  (32) 
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where  C  1,2,3  are  constants  that  depend  on  the  ionization  potential,  and  perhaps  other  hxed 
parameters.  Using  the  ionization  condition  r)(t)  >  H,  one  obtains  for  the  density  of  ionized 
particles 

/rv(t) 

dp  J  g(x,p,  H)dH  (33) 

where  it  is  assumed  the  momentum  of  unionized  particles  is  negligible.  Carrying  out  the 
integration, 

n(x,t )  =  n0(x)  [l  —  e-^]  (34) 

where  n0(x)  is  the  initial  density  of  unionized  particles.  Upon  differentiation  with  respect 
to  time, 

(IT) 

—  =  w(t )  [n0  -  n(t)]  (35) 

which  is  the  correct  macroscopic  ionization  rate,  under  the  stated  assumptions. 


V.  NONLINEAR  PROPAGATION  IN  A  PLASMA  LENS 

A  laser  pulse  propagating  in  plasma  acquires  a  phase  proportional  to  the  plasma  density. 
Plasma  with  density  variation  imparts  a  spatially  varying  phase,  causing  the  pulse  to  refract. 
Thus  with  appropriate  spatial  structuring  the  plasma  can,  in  principle,  be  made  to  mimic 
any  linear,  solid-state  optical  element.  Plasma-based  optical  elements,  being  already  ionized, 
have  the  advantage  of  higher  damage  thresholds,  allowing  their  use  at  higher  intensities  than 
solid-state  elements.  Furthermore,  plasma  optics  can  be  cheaply  and  rapidly  replaced,  for 
instance,  at  the  rep-rate  of  a  gas  jet  or  capillary  [24,  25],  or  flow  rate  of  a  water  jet  [26]. 

A  plasma  lens,  in  which  the  density  profile  of  the  plasma  increases  quadratically  with 
radius,  can  have  enormous  focusing  power  [11,  13].  For  a  density  profile  of  the  form  ne  = 
n0  +  I'Uq't2,  a  short  plasma  lens  imparts  a  quadratic  phase  analogous  to  the  phase  applied 
by  a  thin  lens.  Specifically,  (f)  =  —ifiA/koiv^r2,  where  A  is  the  lens  thickness,  ko  the  pulse 
wavenumber,  wm  =  (2/nrenQ)1^  describes  the  lens  curvature,  and  re  is  the  classical  electron 
radius.  For  an  incident  pulse  with  spot  size  wq  the  effective  /#  =  (l/8)(w^/ Awo)(kowm)2 
provided  /#  >  (A/2 w0).  As  an  example,  we  take  A  =  2-Ji/k0  =  800  nm,  wm  =  15  /mi, 
w0  =  250  /mi  and  A  =  0.5  mm,  and  find  /#  =  3.1. 

A  thick  plasma  lens  can  be  considered  a  truncated  plasma  waveguide  [11,  27].  For  the 
density  profile  above,  the  plasma  waveguide  supports  a  transverse  Gaussian  mode  of  exp(- 
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1)  field  radius  wm.  If  wq  >  wmf  the  pulse  will  undergo  spot  size  oscillations,  reaching  a 
minimum  of  iomin  «  w0(wm/w0)2  after  a  distance  z  ~  7rfc0w^/4,  with  an  effective  /#  = 
(l/4)(wm/wo)koWm.  For  the  same  parameters  as  above,  we  find  /#  =  1.8.  For  both  thick 
and  thin  plasma  lenses  well-established  formation  techniques  can  be  employed.  For  instance 
through  the  gas  ionization,  plasma  heating,  and  hydrodynamic  expansion  driven  by  ~  100 
ps  Nd:YAG  pulse  focused  onto  a  gas  jet  [24], 

The  plasma  lens  configurations  described  above  are  essentially  truncated  versions  of  the 
plasma  channels  that  have  demonstrated  laser  propagation  over  many  Rayleigh  lengths 
[25,  27],  and  have  been  employed  in  channel-guided  laser  wakeheld  accelerators  [28].  At 
moderate  beam  powers,  the  focusing  in  these  plasma  channel  lenses  is  primarily  determined 
by  the  plasma  density  profile.  At  high  beam  powers,  relativistic  self-focusing  enhances  the 
focusing  effect.  A  thin  uniform  plasma  slab  acts  as  a  focusing  lens  if  the  laser  power  is 
above  the  critical  power  for  relativistic  self-focusing  [10].  However,  the  focusing  strength  is 
dependent  on  radial  and  axial  variations  in  the  pulse  intensity,  which  produces  substantial 
aberrations  that  can  significantly  degrade  the  focusing  quality. 

Here  we  examine  the  use  of  a  plasma  lens  to  focus  the  10  PW  beamlines  under  construction 
at  ELI-NP.  Currently  ELI  plans  to  have  two  optical  paths  ending  with  either  an  /A  «  3 
or  f#  fa  20  parabolic  mirror.  A  plasma  lens  placed  within  the  target  chamber  would  allow 
added  flexibility  in  the  focal  geometry.  Additionally,  the  plasma  lens  could  serve  as  a  spatial 
filter  and  help  counteract  any  main  pulse  expansion  resulting  from  pre-pulse  effects.  While  a 
plasma  lens  has  a  higher  damage  threshold  than  a  solid-state  lens,  at  the  extreme  intensities 
of  ELI,  plasma-based  optics  can  still  acquire  aberrations  from  nonlinear  modifications  to 
the  plasma  density.  As  a  result,  neither  the  simple  estimates  for  the  /#  provided  above 
nor  a  weakly  relativistic  approach  will  suffice  [11,  12].  Said  differently,  optimizing  a  plasma 
lens  for  ELI  requires  models  that  can  capture  highly  nonlinear  modifications  to  the  plasma 
density. 

We  take  a  hierarchical  approach  to  optimizing  the  plasma  lens,  using  a  combination  of 
three  models  with  a  varying  degree  of  approximation:  a  computationally  rapid,  nonlinear 
thin  lens  model  based  on  the  beam  propagation  method  (BPM),  ponderomotive  guiding 
center  (PGC)  simulations  based  on  the  modified  paraxial  wave  equation  [29],  and  fully 
electromagnetic  3D  particle-in-cell  (PIC)  simulations  [30] .  Our  starting  point  is  the  thin  lens 
BPM  model,  which  allows  rapid  parameter  scans.  This  model  uses  the  fully  nonlinear  plasma 
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density  perturbation  in  the  adiabatic  limit:  tfwhm  >  1;  where  ujp  =  (47T e2n0/me)1/2  is  the 
plasma  frequency  and  tfwhm  is  the  full-width  at  half-maximum  (FWHM)  duration  of  the 
laser  pulse  [31].  If  the  plasma  lens  is  thin,  f#  >  (A/2w0),  the  laser  pulse  acquires  a  phase 


A 

2/cq7  [ 


klo  +  —rr2  +  Vi7 


wz 


(36) 


on  passing  through  the  lens,  where  7  =  (1  +  | aj_  | 2/ 2) 1//2  and  aj_  =  eA±/mec2  is  the  nor¬ 
malized,  transverse  vector  potential  of  the  pulse.  The  last  term  on  the  right  hand  side  of 
Eq.  (36)  represents  the  transverse  ponderomotive  expulsion  of  electrons  from  the  laser  pulse 
path.  The  q^1  coefficient  includes  the  nonlinearity  responsible  for  relativistic  self-focusing 
in  the  weakly  nonlinear  limit,  |aj_|  <  1.  After  applying  this  phase,  we  can  find  the  pulse 
profile  at  any  distance,  modified  by  the  nonlinear  aberrations,  using  the  beam  propagation 
method. 

The  most  important  parameter  for  the  phenomena  discussed  in  the  previous  sections  is 
the  peak  pulse  intensity.  We  consider  the  ELI  laser  system  and  focal  geometry  with  the 
following  parameters:  A  =  800  nm,  Tfwhm  =  30  fs,  pulse  energy  U  =  200  J,  /#  ~  20, 
and  final  parabolic  mirror  diameter  D  =  0.5  m.  By  itself  this  system  produces  an  intensity 
of  4  x  1021  W/cm2.  As  we  will  see  the  plasma  lens  can  focus  to  intensities  far  surpassing 
this.  The  remaining  parameters  for  optimizing  the  peak  intensity  are  the  thin  plasma  lens 
linear  focal  length  /  =  (l/4)(fc0wm)2(w^/A),  the  location  of  the  plasma  lens  upstream  from 
the  unassisted  focus,  the  plasma  lens  width  for  which  we  choose  A  =  0.5  mm,  and  the 
background  plasma  density  n0.  The  plasma  lens  focusing  power  is  independent  of  n0,  while 
the  nonlinear  phase  associated  with  effects  such  as  relativistic  self-focusing  is  proportional 
to  no.  Consequently  nonlinear  aberrations  can  be  mitigated  without  sacrificing  focusing 
power  by  choosing  the  background  density  as  small  as  possible.  Here  we  use  no  =  1018 
cm-3.  Our  initial  condition  for  the  BPM  model  is  a  pulse  incident  on  the  plasma  lens  with 
the  appropriate  phase  front  curvature  and  amplitude  acquired  by  the  aforementioned  ELI 
parabolic  mirror. 

Fig.  11  displays  the  peak  intensity  as  a  function  of  plasma  lens  location,  vertical  axis, 
and  plasma  lens  linear  focal  length,  horizontal  axis,  for  linear,  7  — >  1,  and  nonlinear  lenses. 
In  both  the  linear  and  nonlinear  cases,  the  peak  intensity  can  be  increased  by  moving  the 
plasma  lens  backwards  (equivalent  to  making  the  lens  larger)  or  by  decreasing  the  focal 
length,  consistent  with  linear  optics.  Even  with  the  nonlinear  aberrations,  the  plasma  lens 
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FIG.  11:  Peak  intensity  as  a  function  of  plasma  lens  location,  vertical  axis,  and  plasma  lens  linear 
focal  length,  horizontal  axis,  for  linear,  7  — >  1,  and  nonlinear  lenses.  Even  with  the  nonlinear 
aberrations,  the  peak  intensity  achieved  by  the  plasma  lens,  ~  2  x  1022  W/cm2,  surpasses  that  of 
the  unassisted  mirror,  ~  4  x  1021  W/crn2. 


FIG.  12:  (a)  On-axis  intensity  as  a  function  of  distance  in  the  speed  of  light  frame,  horizontal 
axis,  and  propagation  distance,  vertical  axis,  (b)  The  transverse  intensity  profile  of  the  pulse  at  a 
propagation  distance  near  where  the  pulse  achieves  its  maximum  intensity. 

achieves  a  peak  intensity  ~  5  times  greater  than  the  unassisted  focus. 

To  examine  the  focusing  more  closely,  we  performed  a  PGC  simulation  of  the  plasma  lens 
focusing  [29].  While  these  simulations  take  longer  than  the  BPM,  they  capture  the  dynamic 
response  of  the  plasma  lens  and  the  non-zero  lens  thickness.  The  same  initial  conditions 
described  above  were  simulated  with  a  plasma  lens  of  focal  length  0.45  mm  located  2  mm 
upstream  from  the  unassisted  focus.  The  laser  pulse  was  initialized  with  a  sin2  temporal 
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profile.  Even  at  this  distance  back  from  the  unassisted  focus,  the  normalized  vector  potential 
is  already  a  sizable  |a_|_|peak  =  8.7. 

Fig.  12(a)  displays  the  on-axis  intensity  as  a  function  of  distance  in  the  speed  of  light 
frame,  horizontal  axis,  and  propagation  distance,  vertical  axis.  The  pulse  reaches  a  peak 
intensity  of  1.3  x  1022  W/cm2:  an  intensity  sufficient  to  ionize  hydrogen-like  argon  or  titanium 
and  «  4  times  greater  than  the  unassisted  focal  intensity.  The  multiple  bright  spots  in 
Fig.  12(a)  result  from  time  slices  within  the  pulse  focusing  at  different  axial  locations  owing 
to  the  varying  degree  of  nonlinearity  encountered.  The  center  of  the  pulse  has  higher  power 
than  the  front  of  the  pulse,  acquires  a  larger  nonlinear  modification  to  its  phase,  and  focuses 
early.  The  back  of  the  pulse,  on  the  other  hand,  focuses  even  earlier  having  encountered 
the  large  ponderomotive  electron  density  modification  driven  by  the  front  of  the  pulse. 
Additionally,  at  a  particular  time  slice,  each  transverse  location  converges  to  the  optical 
axis  at  a  different  rate  producing  the  swath-like  features  trailing  the  bright  spots.  Fig.  12(b) 
shows  the  transverse  intensity  profile  of  the  pulse  at  a  propagation  distance  near  where  the 
pulse  achieves  its  maximum  intensity.  The  brightest  central  spot  corresponds  to  the  center 
of  the  pulse  at  its  peak  intensity.  The  dimmer  spots  to  the  right  and  left  correspond  to  the 
front  and  back  of  the  pulse  before  and  after  their  focuses  respectively. 

While  these  results  are  promising,  the  plasma  lens  focusing  can  be  greatly  improved.  Our 
continuing  research  will  focus  on  further  optimization  by  modifying  the  plasma  lens  to  correct 
for  spherical  aberrations  resulting  from  the  nonlinear  laser-lens  interaction,  structuring  the 
longitudinal  profile  of  the  lens,  using  leaky  lenses  for  mode  cleaning,  and  conducting  full- 
format,  3D  particle  in  cell  simulations  to  examine  non-cylindrically  symmetry  aberrations. 

VI.  RELATIVISTIC  IONIZATION  THEORIES 

Studies  of  ionization  of  an  atom  in  an  electric  held  have  a  long  history  starting  with 
experiments  and  analyses  for  the  case  of  a  constant  electric  held  [32,  33].  A  significant 
advance  was  made  by  Keldysh  when  he  analyzed  the  case  of  a  time-varying  electric  held 
[34],  Since  then  numerous  analytical  and  numerical  approaches  have  been  employed  with 
special  emphasis  on  laser  photoionization. 

Besides  interest  in  photoionization  as  a  fundamental  physical  process  there  are  many 
applications  for  photoelectrons.  Knowledge  of  the  electron  properties,  e.g.,  energy  and 
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momentum  distribution,  is  critical  in  some  of  these  applications  [4,  5]. 

Relativistic  effects  become  significant  either  with  increasing  intensity  or  atomic  number. 
An  example  of  this  is  vacuum  polarization  whereby  electron-positron  pairs  are  created  when 
the  electric  field  approaches  the  Schwinger  field.  Photoionization  of  inner-shell  electrons  in 
high-Z  atoms  is  another  example  where  relativistic  effects  are  important. 

Two  analytical  approaches  are  prevalent  in  studies  of  photoionization.  There  are  several 
interactions  that  must  be  incorporated  in  the  analyses.  The  Coulomb  interaction  with  the 
nucleus  and  the  transition-inducing  interaction  with  the  electromagnetic  field  predominate. 
The  imaginary  time  method  (1TM)  employs  a  Feynman  propagator  to  evolve  the  final  state 
from  the  initial  bound  state.  The  propagator  is  expressed  as  etS^h:  where  S  —  f  dtL  is  the 
action  and  L  is  the  Lagrangian  of  an  electron  in  a  plane  electromagnetic  wave  [35].  The 
action  is  written  as  a  definite  integral  of  L  over  time  and  the  photoionization  amplitude  is 
given  by  minimizing  the  imaginary  part  of  S  using  the  classical  orbits  in  the  barrier  region. 
In  the  second  approach  an  exact  S'-matrix  is  defined  as  an  overlap  integral  between  a  bound 
state  in  the  distant  past  and  a  final  state  in  the  distant  future  [14,  36,  37].  Using  the 
Klcin-Gordon  or  the  Dirac  equations  the  overlap  integrals  are  then  re-written  in  the  more 
familiar  transition  probability  form  involving  the  interaction  inducing  perturbation.  The 
Volkov  solution  (electron  wavefunction  in  the  presence  of  an  electromagnetic  wave)  plays  a 
special  role  in  the  S'-matrix  approach. 

A.  Critique  of  Analytical  Approaches 

The  analytical  approaches  for  determination  of  ionization  rate  incorporate  many  approx¬ 
imations.  Some  of  these  are  listed  here. 

1.  Quasi-  classicality 

The  quasi-classical  nature  of  the  main  exponential  factor,  which  is  common  to  both 
approaches,  requires  that  the  barrier  width  measured  along  the  direction  of  the  electric 
held  (i.e.,  the  transverse  axis  of  the  figure-8  orbit),  be  large  compared  with  the  bound-state 
radius. 


27 


2.  Gauge 


The  ITM  employs  the  radiation  gauge  throughout.  In  analyses  employing  the  S-matrix 
approach  various  gauges  have  been  considered.  However,  different  gauges  lead  to  different 
results.  The  rationale  for  the  choice  of  gauge  has  been  discussed  in  Klaiber.  It  is  known  that 
ionization  rates  based  on  ITM  are  in  good  agreement  with  experimental  results.  Therefore 
it  is  argued  that  the  appropriate  gauge  for  the  ^-matrix  approach  is  the  one  that  recov¬ 
ers  the  ITM  rates.  This  is  the  length  gauge  in  non-relativistic  theory  which  generalizes 
to  the  Goppert-Mayer  gauge  in  relativistic  theory.  There  is  still  the  issue  of  partitioning 
the  Hamiltonian  operator  into  the  bound-state  Hamiltonian  and  the  interaction  Hamilto¬ 
nian.  The  Coulomb-corrected  dressed  ionization  rate  corresponds  to  the  choice  where  the 
bound-state  Hamiltonian  partially  includes  the  electromagnetic  field  while  the  interaction 
Hamiltonian  includes  the  remainder.  The  former  takes  into  account  the  interaction  of  elec¬ 
tron  spin  with  the  electromagnetic  field  due  to  i)  bound-state  energy  level  shifts  (i.e. ,  Zeeman 
splitting)  and  ii)  electron  spin  precession  in  the  laser  field,  while  transition  to  the  final  state 
is  clue  to  the  latter. 


3.  Nuclear  Coulomb  Tail 

Except  for  negative  ions  (short-range  potential)  the  nuclear  charge  has  a  significant  effect 
on  the  wavefunction  of  the  ejected  electron  and  therefore  on  the  ionization  rate  and  is 
incorporated  in  both  methods  through  a  multiplicative  factor.  The  approximations  involved 
in  analytical  S'-matrix  approaches  limit  their  validity  to  ionization  potentials  <C  0.5  MV. 

B.  Ab  Initio  Simulations 

In  light  of  the  limitations  discussed  above,  an  important  part  of  the  overall  effort  to  de¬ 
scribe  a  process  such  as  xLIPA  lies  in  benchmarking  the  ionization  rate  laws  used  in  the 
particle  tracking  calculation.  For  this  purpose  NRL  has  developed  a  suite  of  ab  initio  simu¬ 
lation  models  which  solve  various  quantum  mechanical  wave  equations  for  a  charged  particle 
in  an  arbitrary  electromagnetic  field.  The  wave  equations  that  have  been  incorporated  into 
the  suite  are  the  Schrodinger,  Pauli,  Klein-Gordon,  and  Dirac  equations.  The  entire  suite 
of  models  takes  advantage  of  a  combined  OpenCL/MPI  programming  model  that  scales  to 
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large  numbers  of  GPGPU  devices  operating  in  parallel.  A  subset  of  the  models  and  results 
are  documented  in  the  literature  [38-40]. 


VII.  CONCLUSIONS 

The  xLIPA  scheme  of  free  space  electron  acceleration  promises  to  produce  highly  direc¬ 
tional,  GeV  class  electrons,  by  means  of  direct  acceleration  in  the  fields  of  a  multi-PW  class 
laser  pulse.  In  the  process,  nuclei  of  moderate-Z  atoms  such  as  argon  or  titanium  are  pro¬ 
duced  through  tunneling  ionization.  The  energy  of  the  electrons,  and  charge  of  the  nuclei, 
would  be  record-setting  for  free  space  acceleration  and  optical  tunneling  ionization,  respec¬ 
tively.  In  the  event  that  multi-EW  laser  pulses  should  become  available,  radiation  reaction 
effects  lead  to  observable  effects  in  the  momentum  distribution.  One  way  to  reduce  the  laser 
power  necessary  for  this  is  to  improve  the  laser  focus. 

In  order  to  improve  the  flexibility  and  efficacy  of  potential  experimental  configurations, 
a  plasma  lens  is  proposed  as  a  combined  final  focusing  element  and  target.  In  the  lin¬ 
ear  approximation,  the  enhancement  in  irradiance  due  to  the  plasma  lens  may  be  orders 
of  magnitude.  When  nonlinear  effects  are  taken  into  account,  lens  aberrations  limit  the 
enhancement  to  factors  of  several.  Investigations  are  underway  to  determine  whether  the 
aberrations  can  be  corrected  by  tailoring  the  plasma  lens  density  profile.  In  any  case,  the 
lens  density  profile  weights  the  initial  particle  positions  in  favor  of  high  final  energy. 

The  xLIPA  scheme  has  several  benefits,  both  in  terms  of  applications,  and  fundamental 
physics  interest.  One  possible  application  is  as  an  injector  for  a  staged  laser  wakefield 
acceleration  system.  Another  is  as  a  source  of  high  energy  electrons  to  be  used  in  driving  an 
x-ray  free  electron  laser,  or  other  advanced  light  source.  The  fundamental  physics  of  xLIPA 
is  rich,  consisting  of  relativistic  tunneling  ionization,  charged  particle  dynamics  in  extreme 
fields,  and  eventually  radiation  reaction. 
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